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Abstract 

One of the most fundamental concepts of evolutionary dynamics is the "fixation" 
probability, i.e. the probability that a mutant spreads through the whole pop- 
ulation. Most natural communities are geographically structured into habitats 
exchanging individuals among each other and can be modeled by an evolution- 
ary graph (EG), where directed links weight the probability for the offspring of 
one individual to replace another individual in the community. Very few exact 
analytical results are known for EGs. We show here how by using the tech- 
niques of the fixed point of Probability Generating Function, we can uncover a 
large class of of graphs, which we term bithermal, for which the exact fixation 
probability can be simply computed. 

Keywords: Evolutionary graphs, fixation probability, fitness, probability 
generating functions. 



1. Introduction. 

Evolutionary dynamics is a stochastic process due to competition between 
deterministic selection pressure and stochastic events due to random sampling 
from one generation to the other. One of the most fundamental concepts of 
evolutionary dynamics is the fixation probability, i.e. the probability that a 
mutant spreads and takes over the whole community (Patwa and Wahl |12|). In 
the framework of the Moran model (Moran [ll|) for a well mixed population, 
where an individual's offspring can replace any other one in the community, the 
fixation probability is 

1 _ r -m 

*f = t^-m (i) 

where M is the size of the community, mo the original number of mutants and 
r the relative fitness of the mutants. A similar result was reached by Kimura 
(Kimura |7|) for the Fisher- Wright model under the diffusion approximation. 

The idea of well mixed population is however far from realistic. Natural com- 
munities are geographically extended and subdivided into patches that exchange 
individuals (figure QJi) . Maruyama (Maruyama @, was the first to cast the 
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Figure 1: Evolutionary Graphs, (a) individuals are spread in space. A node i can send its 
offspring only to connected nodes j with a probability rriij. A non-structure population can 
be considered as a fully connected graph with uniform migration probability rriij = l/N. (b) 
The star configuration, which can be proved to differ from a non-structured population for its 
fixation probability. 

problem of evolutionary dynamics into a Moran process on graphs (or islands, 
in his terminology) and under harsh approximations, concluded that the fixa- 
tion probability does not depend on the population structure. The first formal 
proof that Maruyama's results are not always correct was given by Lieberman 
et aZ(Lieberman et al. [f|) who showed that for a Moran process on a star graph 
(Figure [TJ)), the effective fitness of the mutant, in the large population limit, 
can be enhanced to r 2 due to topological effects. In order to do so, Lieberman 
et al. considered evolutionary graphs (EG) where nodes contain exactly one 
individual, whether wild type or mutant, connected by directed links represent- 
ing the geographical (or social) connectivity. Lieberman et al. also extended 
their results to the funnel topology with K layers where the effective fitness, 
in the limit of large population, can be amplified to r K , but provided only a 
sketch of the proof. Beyond the special cases considered by Lieberman et al, 
very few exact analytical results are known. A review of the present state of 
known results is given by Shakarian et al(Shakarian et al. fl3|). 

Most of the results of the EG are obtained through Monte Carlo numerical 
simulations. These simulations however scale as 2 M where M is the number 
of nodes. A new numerical scheme has been proposed (Barbosa et al. [2j) to 
accelerate the speed of these simulations, but the computation of the fixation 
probability of large graphs is still very time consuming. 

It would therefore be important to know the exact fixation probability of a 
large class of graphs that can be used as an approximation of closely related 
graphs or as a benchmark for assessing the progress in numerical simulation 
schemes. This is the aim of the present article. 

We recently proposed a new method (Houchmandzadeh and Vallade Q), 
based on the fixed points of the time dependent Probability Generating Function 
(fp-PGF) which can efficiently approximate the fixation probability of large, 
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arbitrary graphs by solving only a system of M second order algebraic equations. 
We show in the present article that the fp-PGF method can also be used to 
derive exact results for a large class of graphs that we call bithermal, which are 
an extension of the isothermal graphs considered by Lieberman et al. 

The temperature of a node is related to the imbalance between the sum 
of the weights of the incoming and outgoing links. In isothermal graphs, all 
nodes are balanced and have the same temperature T = 1. Bithermal graphs 
are bipartite, with two kinds of nodes at either temperature Ta or Tb- The 
star topology (figure [TJd) is one particular example of such graphs, some other 
particular examples are shown in Figure [5] . We show here that the fixation 
probability of these graphs is a simple rational function of the fitness r and of 
the number of nodes Ma and Mb in each class 

K f = f(M A ,M B ,r) 

The exact expression of this function is given by equation (|2ip and its plot by 
figure When there is the same number of nodes in each class (Ma = M B ), 
bithermal graphs become isothermal and the function / above is equal to the 
Moran expression (fl}. When the imbalance between the number of nodes in 
each class is large (Ma S> Mb or Ma <C Mb), the fixation probability tends 
toward or 1 — 1 /r 2 , depending on the nature of the Moran process (birth-death 
or death-birth). 

This article is organized as follow. In the next section, we recall the con- 
tinuous time stochastic process of Moran on graph and its associated Master 
equation ; the third section is devoted to the Probability Generating Function 
method ; In the fourth section, we apply these results to the bithermal graphs 
and give their exact fixation probability. The last section is devoted to some 
generalizations and conclusions. 

2. Continuous time Moran process on graph. 

Consider a community of M individuals, which can either be wild type (WT) 
with fitness 1 or mutant with relative fitness r = 1+s. The individuals are spread 
spatially and the progeny of an individual i can replace individual j according 
to a connectivity map. The connectivity map can be envisioned as a graph, 
where each node i contains exactly one individual, either mutant or wild type ; 
the weight of a link specifies the probability for the progeny of an individual 
at node i to replace an individual at node j (Figure |TJl) . The coefficients m^- 
are collected into a connectivity matrix m. As the number of individual is 
fixed, it is sufficient to specify the number of mutants (0 or 1) on each node 
n = (m, ri2, ■■■Tim) at a given time to have complete information about the 
system at this time. We consider here a continuous time model where birth 
(or death) events occur randomly with rate \i (Houchmandzadeh and Vallade 
[5]). The probability density for a node i to decrease or increase its number of 
mutants by one unit during a time interval dt is (Houchmandzadeh and Vallade 
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Figure 2: Some examples of bithermal graphs; for simplicity, each link represents two directed 
links. A bithermal graph is constituted of two kinds of nodes at respective temperature Ta 
and Tg , indicated here by their colors. Each A node is connected to a subset of B nodes and 
vice versa, (a) a (4,12,6,2) generalized star, consisting of 4 central nodes and 12 peripheral 
nodes, each central node connected to 6 peripheral and each peripheral connected to two 
central nodes ; (b) an example where the combined effect of the number of links and their 
weight makes the graph bithermal ; (c) A symmetric 2 dimensional, bithermal crystal with 
periodic boundaries; 



4 



1) 

Wr(m) = (^/M)ni^m fci (l - n k ) (2) 



M 



k=l 

M 

Wf{m) = {n/M)r(l-ni)J2mkin k (3) 

fe=i 

It is crucial at this step to distinguish between two kinds of Moran processes 
(Antal et al. In the first case which we call D-B (Death first and then 

replacement, also called Voter Model), a death occurs first on a node, then 
immediately one connected node duplicates and send its progeny to this node. 
Equation @ is therefore the probability density that a mutant dies at node 
i during dt, and is replaced by the progeny of a WT on a connected node k. 
Equation is the probability density that a WT dies at node i and is replaced 
by a mutant on a connected node. Without loss of generality, the mutant's 
advantage r is included in this line, either as a decreased mortality or a better 
replacement success once a death event has occurred. 

In the case of B-D processes (also called Invasion Process), a birth occurs 
first on a node fc and the progeny is then sent to a connected node i to replace 
the local resident. The transition probabilities are still expressed by the same 
equations (|2I3|) . but the quantity /i now denotes the birth rate. 

Although the rate equations (|2l3p are similar for these two processes, the 
normalization conditions of m ki coefficients are different : 

M M 

J2 m ik = 1 (D-B); ^ m kl = I (B-D) Vfc (4) 

i=l i=l 

The temperature of a node is defined for these processes as 

M M 

m kl = T k (D-B) ; £ m lk = T k (B-D) Vfc (5) 

i=l i=l 

Because of the normalization constraint we must have Y]j—i Tj = M for 
both processes. This means that if some nodes are cold (T < 1), others must 
be hot (T > 1). 

The Moran process is a one-step stochastic one, where during an infinites- 
imal interval dt, only one birth or death event can occur. Equations (|2I3|) 
are transition probabilities between states n on the one hand and states a^n = 
(ni, ...n,- — 1, ..Mm) and ajn = (m, ...n,- + l, ...njw) on the other. The probability 
P(n, t) of observing state n at time t obeys the Master equation 

^P< ^ , - = W(m -> n)P(m, t) - W(n -> m)P(n, t) (6) 

r i 
{m} 

The EG process we are considering has two absorbing states 1 = (1,1, ...1) 
and = (0, 0,...0). Once a mutant has invaded all the nodes or has been 
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eliminated from all of them, it is fixed or lost and there is no further evolu- 
tion (until a new mutant appears by random mutation) : W(l — > m) = and 
W(0 —> m) = as can be deduced from equations (|2I3I) . Note that the probabil- 
ity of reaching state 1 from an initial state n, i.e. the fixation probability 7r(n) 
can in principle be found using the Kolmogorov's backward equation (Ewens 

1): 

(7r(n) - 7r(m)) W(n ->■ m) = (7) 

m 

tt(0) = ; tt(1) = 1 (8) 

which is a set of linear equations in the unknowns 7r(n). The direct resolution of 
the above set of equation however can be attempted only in special cases. The 
best example of a direct solution is the unstructured population, where the graph 
is fully connected and my = 1/M ; the EG dynamics can then be mapped into 
a biased one-dimensional Brownian motion and solved by standard techniques 
(Ewens 0]), which yields the well known result ([1]) . Other cases, such as the 
star topology in the limit of large population considered by Lieberman et al 
(Lieberman et al. [8j|), use such careful mapping. The mapping method however 
is hard to generalize. 



3. The fp-PGF method. 

Computation of the fixation probability can be simplified if instead of the 
linear system ([?])■ we use the dynamics of the Probability Generating Function 

0(z,<) =^P(n,f)z n 

where the variable z = (z\, Z2, ...zm) is conjugate to n = (n\, ri2, tim) and 
z n = z™ 1 z 7 2 12 ...z^f 1 . Time is measured in generation time units M/fj,. Note that 
</>(0,t) = and 0(1, f) = 1. From the Master Equation ©, we can derive the 
dynamics of the PGF (Houchmandzadeh and Vallade [6]) which reads: 

d( t ) t ( \ d $ \^ < \ d2 <t> /r^ 

^ = £ /fe(z te~?/ i ' fe(z) ^ () 

k—l i.k—1 



where for D-B processes, 

/fe(z) = z k \ y ]m kl (zi - 1)1 - r~ l (z k - 1) (10) 
Si,fc(z) = m ki (z l -l)(z l -r^ 1 )z k (11) 



( M \ 

I y^kijzj - 1) J - r l {z k -l) 



For the B-D process, the first order term is slightly different and reads 

M 

fk(z) = z k I y^ m kl (zi -1) 1 - (T k /r)(z k - 1) 



^m kl (zi - 1) J - (T k 



G 



Solving the PGF equation © would seem at least as formidable as solving 
directly the Master Equation ([6]). However, for the computation of the fixation 
probabilities, we are only interested in the large time limit t — >• oo. At large 
time, the mutant is either fixed or lost, therefore the stationary solution (f> s (z) 
to which the PGF converges is simply 



,(z)=7r +7r/JJzj (12) 



i=l 



where ttq and ttj are the loss and fixation probabilities and implicit functions 
of the initial conditions. It can also be checked manually that (fT2|) is indeed a 
solution of ([9]). The problem of finding tto and 717 becomes trivial if the PGF 
possesses a fixed point £ = (£1, ...£a/) such that 



0t 

In this case, we have 



= 



A I 



<KC,0) =0(C,OO) = TTO + 7T/ IJ C< 

i=l 

and, as 717 + 7To = 1, 

i-0(C,Q) 
7T/ = : — =m — 



1 - KU c 

As the quantity 0) is known from the initial conditions, finding a fixed point 
of the PGF determines entirely the fixation probability. Note that once a fixed 
point has been found, the fixation probability for any initial condition can be 
trivially computed. 

For the initial condition of the mutant appearing at random with probability 
1/M on one node, <j>(zi, ...ZM',t = 0) = (1/M) ^ Zi and therefore the fixation 
probability is given by 

*' = i-n,c, (13) 

The condition for a point ( to be a fixed point is 

A(C) = ovfc (w) 

9i,k{Q+9kAQ = 0Vj,fc (15) 

Whether such a fixed point exists or not depends on the connectivity matrix 
rriik- For an isothermal graph where all the nodes k have temperature Tk = 1, 
it is easy to check that £ = r _1 l is a fixed point : 

fkir- 1 !) = r- 1 (r- 1 -l)-r- 1 (r- 1 -l)=0 
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M A M B 




Figure 3: The connectivity matrix of a 2-thermal graph for a D-B process, where nodes of 
type A are numbered from 1 to Ma and nodes of type B from Ma + 1 to Ma + Mg . For a 
B-D process, the sum of the elements of a row is equal to unity, and the sum of the elements 
of a column is equal to Ta or Tg. 

and the fixation probability (I13[) of the isothermal graph is equal to the result 
(fTj) for unstructured populations. We see here how easily this theorem can be 
obtained from the fixed point of the PGF. 



4. bithermal graphs. 

We now consider a subset of bipartite graphs that we call bithermal 
for which the fixation probability can be determined in algebraic closed form. 
In these graphs, there are Ma nodes of type A at temperature Ta and Mb 
nodes of type B at temperature Tg. More over, we require that two nodes 
can be connected only if they are at different temperatures, and if a node i is 
connected to a node j, then j is also connected to i. Of course, we suppose that 
the graph is connected, i.e. there is always a path from any node i to any node 
j. With appropriate numbering of the nodes, the connectivity of such graph is 
a block matrix of the form 



(figure [3]). The star topology (figure [TJd) is such a graph with Ma = 1 central 
and Mb peripheral nodes. For a D-B process, all the elements of the star's /3 
block are equal to 1/Mb and all the elements of a block are equal to 1. Therefore 
Ta = Mb and Tb = 1 /Mb- The star graph can be generalized to the case where 
the number of central nodes Ma > 1 (figure [2^) . For a general bithermal graph 
the number of B nodes to which an A node connects and the weights of these 
links can be arbitrary, as long as the constraints on temperatures are respected. 
Note that Ta and Tb are not independent. By summing over all the elements 
rriki of the j3 or a block of the connectivity matrix m, we have for a D-B process 

M A = M B T B ; M B = M A T A 
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for a B-D process, the role of Ma and Mb are exchanged, but in both cases, we 
have 

T A T B = 1 (16) 

We now search for the bithermal graphs which have an exact fixed point. 
Following the example of the isothermal graphs, we look for a fixed point £ — 
(Ci,— Cm) where Q = (a if i £ A and Q = Cb if i € B . 

Let us first consider the case of D-B processes. Using condition (fl4|) for 
k G A or B we obtain a set of two algebraic equations : 

MO = a(Cs-i)T A -r- i (a-i) = o 

MC) - Cs(a-l)T B -r- 1 (C s -l) = 

the solution of which is given by 

(a = u(r, Ma/M b ) (17) 

Cs = u(r,M B /M A ) (18) 
l/r + T 

<r,T) = (19) 

This solution also satisfies condition (fT5j) if for fc g A and i £ B, we have the 
following relation between the coefficient of the connectivity matrix m : 

m u (Cb - 1)(Cb - t^Ka 



m lk (Ca ~ 1){(a - r-^CB 

1 _ M A 
T a M b 

which we can express as a relation between the two blocks j3 and a of the 
connectivity matrix m: 

a = T A p J (20) 

We observe here that this is the sufficient condition to form exactly solvable 
bithermal graphs : form an Mb x Ma matrix /3 where the sum of elements 
in each column is 1 and the sum of elements in each row is 1/Ta — Ma/ Mb; 
form the block matrix m from (3 and a = Ta/3 1 . The f3 block contains MaMb 
coefficients subject to Ma + Mb summation rules, so the number of bithermal 
graphs with exact solutions is indeed very large when Ma or Mb are large. 

An important subset of bithermal graphs that always has an exact fixed 
point is a set we call symmetric bithermal graphs. For these graphs, all the 
existing links from a node A to a node B (respect. B to A ) have the same 
weight itiab (resp. tuba)- The generalized star graph (fig[5Ji) belongs to this 
subset. For these sets, the weight of a link is determined only by the number 
of connected nodes, and each A (resp. B) node has always the same number of 
neighbors. The symmetric subset automatically satisfies all the constraints and 
always has exact fixed points. 
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Once the fixed point is known, the fixation probability is easily determined 
from equation (|13l) : 

_ 1 - (MaCa + M b Cb)/{M a + M B ) , , 

77 f i ^Maz-Mb ^ > 

For a B-D process, the computation of the fixed point follows exactly the 
same line of argument. The result is obtained by permuting Ma and Mb '■ 

Cf D = u(r, Mb/Ma) ; Cf D = u(r,M A /M B ) (22) 

The fixation probability is given by the same expression (|2"Tj) . Note that D-B 
and B-D processes act in opposite directions. Compared to a non-structured 
population of the same size, a bithermal D-B process acts as a suppressor of 
selection, lowering the fixation probability where the B-D process is an amplifier 
of selection, increasing the fixation probability. The maximum 7T/ for a D-B is 
obtained for Ta — 1 and is equal to expression ([T]) ; the maximum it/ for a B-D 
process is obtained for Ta — > oo and is equal to 1 — r~ 2 . A plot of both fixation 
probabilities and their comparison to numerical simulations is shown in figure 

m. 

We stress that for bithermal graphs, the details of the connectivity are not 
important : the fixation probability depends only on the total number of A and 
B nodes. Consider for example the symmetric generalized star (Ma, M B ,p, q) 
where each A node is connected to p nodes of type B and each B nodes to q 
nodes of type A. For a fully connected generalized star, p = Mb and q = Ma 
; an example of partially connected generalized star is given in figure [2^. We 
emphasize that for generalized stars, the fixation probability does not depend 
on the detail of the connectivity p, q, a result which would be hard to predict 
by other methods. We also note from numerical simulations that the fixation 
time is also only a function of Ma and Mb and does not depend on the detail 
of the connectivity. 

Finally, we note that even when the constraint (|2U)) is not respected and 
a =/= TaP 1 , the fixation probability of bithermal graphs is well approximated 
by expression (|21l) . In this case, the point £ computed from equations (| 171181) 
is only a quasi-fixed point, but as we have shown earlier(Houchmandzadeh and 
Vallade [6]), for large communities, quasi- fixed points can be used for a good 
approximation of the fixation probability. It can be observed in figure (|3Jd) that 
the numerical errors of fixation probabilities of connectivity matrices having 
exact fixed points (left of vertical line) or only quasi-fixed point (right of vertical 
line) are of the same magnitude, for a system as small as M = 100. We also 
note from numerical simulations that fixation time, is only 



5. Discussion and Conclusion. 

The approach we presented above can be generalized in various directions. 
We have restricted our approach to the case where each node contains only one 
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Figure 4: Fixation probability izf of bithermal EGs of size M = 100 and r = 1.25. (a) TTf 
as a function of the number Ma of nodes of type A for a generalized {Ma, Mb) star, circles 
(black) and triangles (red) correspond to numerical resolution of the fixation probabilities for 
B-D and D-B processes. Solid lines correspond to exact solutions given by equation (ST) . The 
dashed lines correspond to 7r/ = l — 1/r and irj = 1 — 1/r 2 . Numerical simulations were 
performed by using a Gillespie Algorithm(Gillespie 4]) ; for each point, 10 5 stochastic paths 
were generated and used to compute the fixation probability. (b)Fixation probability 7Tf for 
random bithermal EG for different values of Ma (circles 5, squares 10, X 15, diamond 20, 
triangle up 30, triangle down 45) and D-B processes. To the left of the vertical line (#index 
<50) : Each point corresponds to a random bithermal connectivity matrix where the two 
blocks are related through a = TaP 1 ■ To the right of the vertical line (#index >50) : For 
each point, the two blocks of the bithermal connectivity matrix are random and a ^ Ta0 t . 
Solid lines indicate the theoretical values computed by equation (12 H . 
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Figure 5: Generalization to bi-level graphs, where each node can be considered as composed 
of N individuals, conserving the same topology, i.e. the weight of new links is downscaled by 
a factor N : m£ • = rtiij/N 

individual. This restriction can be relaxed and we can let each node contain 
a number of individuals N > 1 (figure [S]) . This is equivalent to the original 
island model of Maruyama(Maruyama [9|) or can alternatively be envisioned as 
a bi-level graph as defined by Shakarian et a/(Shakarian et al. [HI), where each 
node is mapped into N sub-nodes conserving the original topology. As we have 
shown earlier (Houchmandzadeh and Vallade |6j), this parameter does not alter 
the expression of the PGF and the fixed points are still computed by the same 
equations. The fixation probability is slightly modified and reads 

_ 1 - (MaU + M b Cb)/(M a + M B ) 

n f -, /-NM A +NM b 

The approach can be extended even further and allows for different numbers 
of individuals for each node k (Houchmandzadeh and Vallade [(|). 

Another direction toward which this work can be extended is the n— thermal 
graphs. Here we have considered only bithermal graphs composed of two types 
of nodes which we can represent by an A — B topology. In principle, we can 
generalize the method to graphs containing P types of nodes Oi,...,Op : The 
Mj nodes belonging to type Or have the temperature Tj and a hot class can 
only be connected to cold classes and vice versa. We could in principle form 
a polymeric topology such as 0\ — 0%... — Op, branched systems, closed rings 
and so on. The exploration of these topologies implies an analytic study of the 
roots of algebraic equations of order 2P and is beyond the scope of the present 
article. 

To summarize, we have obtained exact analytical results for a wide class 
of graphs that we have called bithermal in the field of Evolutionary Graph 
Theory. EGT is a cornerstone for our understanding of evolution, because 
natural population are always geographically extended and cannot be a priori 
approximated as "well mixed". Exact results in EGT however have been hard to 
obtain because in each case, a mapping into a one-dimensional Brownian motion 
has to be constructed ; whether such a mapping exists or not for a particular 
problem is not a trivial problem. The method we develop is radically different 
: by using the continuous time version of the Moran model and the dynamics 
of the PGF, we reduce the problem of finding an exactly solvable model into 
finding the roots of algebraic equations. We have illustrated the power of this 
dynamical method through our study of bithermal graphs. We believe that the 
method we have presented can be a powerful tool to get exact results for the 
fixation probability of more complex evolutionary graphs. 
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